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CHAOS IN THE LORENZ EQUATIONS: 
A COMPUTER-ASSISTED PROOF 

KONSTANTIN MISCHAIKOW AND MARIAN MROZEK 



Abstract. A new technique for obtaining rigorous results concerning the 
global dynamics of nonlinear systems is described. The technique combines 
abstract existence results based on the Conley index theory with computer- 
assisted computations. As an application of these methods it is proven that 
for an explicit parameter value the Lorenz equations exhibit chaotic dynamics. 



Introduction 

The purpose of this note is to briefly describe a new technique for obtaining rig- 
orous results concerning the global dynamics of nonlinear systems. The technique 
combines abstract existence results based on topological invariants (the Conley 
index) with finite, computer-assisted computations necessary to verify the assump- 
tions of the theorems in a concrete example. There are at least three new aspects 
of this technique. It applies to concrete differential equations; it can provide a 
relatively strong description of global dynamics (in terms of semiconjugacies); and 
the necessary computer-assisted computations are small enough to be performed 
on currently available computers. 

To focus the presentation of the ideas, an outline of a proof that the Lorenz 
equations, 

i ^ s{y-x), 
(1) y = Rx — y — xz , 

z = xy-qz, 

contain chaotic dynamics for a prescribed open set of parameter values will be 
presented. 

Let / : R" —>■ R" be a homeomorphism. For N C R" , the maximal invariant 
set of N is defined by Inv(7V , /) = {x e iV | /"(x) e iVV n e Z} . 
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Theorem 1. Let 

P := {{x , y , z) \ z = 53}. 

For all parameter values in a sufficiently small neighborhood of {s , R, q) = (45 , 54 , 10) ,| 
there exists a Poincare section N c P such that the Poincare map g induced by 
(1) is Lipschitz and well defined. Furthermore, there exists a G N and a contin- 
uous surjection p : liiv{N , g) ^ Tt2 such that 

po g'^ = a o p 

where o" : E2 ^ E2 is the full shift dynamics on two symbols. 

The choice of parameters in Theorem 1 was dictated mainly by an attempt to 
minimize the necessary computations. Preliminary estimates show that analogous 
computations for the classical choice of parameters (.s , R , q) = (10, 28, 8/3), 
though substantially more complex, are also in the range of currently available 
computers (computations in progress, see [4]). However, it must be emphasized 
that Theorem 1 shows only the existence of an unstable invariant set which maps 
onto a horseshoe. Though the set may lie within a strange attractor, we do not yet 
have sufficiently strong abstract results to prove that the whole attractor is chaotic. 

Let us mention some other attempts to reduce the question of chaos in the Lorenz 
equations to a finite computation: an fi-explosion approach [7] and a shooting 
method approach [1]. Contrary to the presented method, in both cases the au- 
thors have not yet performed the actual computations.^ Moreover, as mentioned 
before. Theorem 1 provides a description of the dynamics in a neighborhood of an 
explicitly presented parameter value, i.e. {s , R, q) = (45 , 54 , 10) . This should 
be contrasted with the other methods where the conclusion is that explicit dynam- 
ics occur for some unknown parameter value within a specified range (shooting 
method) or for a sequence of values tending to a given parameter value (the O- 
explosion method). 

The proof of Theorem 1 has five distinct components: 

1. algebraic invariants based on the Conley index theory which guarantee the 
structure of the global dynamics, in this case the semi-conjugacy to the full 

two shift; 

2. an extension of these invariants to multivalued maps; 

3. a theory of finite representable multivalued maps which, when combined 
with the above-mentioned algebraic invariants, serves to bridge the gap 
between the continuous dynamics (in this case the Lorenz equations) and 
the finite dynamics of the computer; 

4. the numerical computations of the finite multivalued map of interest; 

5. the combinatorial computations of the Conley index for the multivalued 
map. 

The rough idea of the general scheme of the proof is as follows. Choose a potential 
isolating neighborhood N for the Poincare map g . Select a finite representable 
multivalued map Q (definitions follow) such that Q is an extension of g , i.e. 
.9(2;) G G{x) for all X G dom Q{x) . Perform a computer calculation both to 
determine Q and to check whether N is an isolating neighborhood for Q . Theorem 
3 guarantees that whenever A'' is an isolating neighborhood for ^ , it is also an 



^ Since completion of the original version of this paper, numerical calculations based on modi- 
fications of the ideas of [1] have been completed [8]. 
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isolating neighborhood for g and the Conlcy indexes coincide. If N is an isolating 
neighborhood, then the Conley index of N under Q is computed (this is again 
a finite computation and can be done by computer). This then determines the 
Conley index for the set under the original Poineare map g and allows one to 
verify the assumptions of Theorem 2, which yields Theorem 1. 

If A'' is not an isolating neighborhood of Q , then one is free to choose a new 
multivalued map which is a better approximation of g and to repeat the compu- 
tations. Theorem 4 implies that if N is an isolating neighborhood for g , then N 
is an isolating neighborhood for a sufficiently refined choice of ^ . In other words, 
if N is an isolating neighborhood under g , then given sufficient computing power, 
this scheme will provide a proof of this fact and, in addition, the values of the 
Conley index of N under g . 

In the sections below each of the five components of the proof will be discussed. 
Most of the theoretical results hold true in greater generality; however, for the sake 
of exposition the simplest setting has been chosen. 



Again, let / : R" — > R" be a homeomorphism. A compact set N is called an 

isolating neighborhood if Inv(A'^ , /) C int A'^ . To determine the Conley index of an 
isolated invariant set for a map, one begins with an index pair {N , L) . The details 
of the definitions which follow can be found in [5]. For the purposes of this paper 
it is sufficient to know that A*" is an isolating neighborhood, L is its exit set, and 
the map / induces a homomorphism called the index map on the cohomology of 
the index pair, i.e. : H*{N , L) H*{N , L). The cohomological Conley index 
of an isolating neighborhood under / is given by 



where CH*{lnv{N , /)) is the graded module obtained by quotienting H*{N , L) 
by the generalized kernel of I J and x*(Iiiv(A'^ , /)) is the induced graded module 
automorphism on CH*{lny{N , /)) . 

Let N = Nq U Ni be an isolating neighborhood under / where A^o and A''i 
are disjoint compact sets. For k , I ~ 0,1, let A^fc; = Nk Ci f{Ni). Let Sk ■= 
Inv(A^fcfc , /) and Sik ■= lnv{Nkk U Nki U Nu , /) . The following result is a special 
case of [3, Theorem 2.3]. 

Theorem 2. Assume that 



and that X*{^ik) is not conjugate to X*{^k) © X*{^i) ■ Then there exists d € N 
and a continuous surjection p : lnv{N , /) ^ S2 such that 



Recall that a multivalued m,ap from R" to itself is a function : R" 'P(R") \ 
{0} from R" to the power set of R" . A continuous function / : R" R" is a 
selector iom if f{x) € T{x) for all x € R" . In the other direction, a multivalued 



1. Isolating neighborhoods and chaos 



Con*(Inv(Ar, /)) = (cH*iIay{N , /)) , x*(Inv(Ar, /))) 




p O = £7 O p 

where cr : S2 ^ S2 is the full shift dynamics on two symbols. 



2. Index theory for multivalued maps 
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function is an extension of a continuous function / : R" R" if /(x) G T{x) 
for all X G R" . Let be a multivalued function on R" . For A C R" let 
^(A) := Ua;GA-^(^) ^^'l define, recursively, := ^(^(^)) . In this way 

T defines a multivalued discrete scmidynamical system on R" . A Conley index 
theory for such systems has been developed in [2], and the basic ideas are as follows. 
Given B c R" , its inverse image is ^~^(B) := {a; e R" | C -B} and its 

weak inverse im.age is T*~^{B) {a; e R" | J^{x) r\ B ^ 0}. Given iV C R" , 
the invariant set of N is given by 

Inv(iV , J^) := {x € N IB-y^ -.Z ^ N such that 7^(0) = , and 

The diameter of ^ over N is the number 

diam atJT := sup{l|z - y\\ \ z , y e Tix)}. 
xeN 

The diameter of over its domain will also be called the size oi . N is called 

an isolating neighborhood under if 

B(Inv(iV , jf) , diamjv C int (N). 

In order to define the Conley index for ^ in , some conditions (admissibility) 
must be met (see [2]). For the sake of this paper it is enough to mention that 
convex-valued maps with continuous selectors are admissible. 

Theorem 3. Let T be an admissible multivalued map and f a selector of T ■ If 
N is an isolating neighborhood for , then it is an isolating neighborhood for f . 
Furthermore, 

Con* (lnv(Ar, /)) f« Con* (lnv(iV, JF)). 

The importance of this theorem is that it implies that if the Poincare map g 
given by the Lorenz equations is replaced by an extension Q (i.e. a multivalued 
map such that g{x) G Q{x)) such that N is an isolating neighborhood for Q , then 
any index information obtained for Q is valid for g . 

The next theorem is a typical convergence theorem and states that any index 
information of g can be determined by a sufficiently small extension Q . 

Theorem 4. Let N be an isolating neighborhood for f : R" R" , a Lipschitz 
continuous function. Let {j^n} be a family of extensions of f such that J^n ~> / • 
Then for n sufficiently large, N is an isolating neighborhood for T ■ 

For the purposes of computation it is convenient to use the following special form 
of an isolating neighborhood. N is an isolating block for if 

(2) B{j^*-\N) n N n J^{N) , diamN J^) C int (iV). 

Notice that contrary to the notion of an isolating neighborhood, it uses only a finite 
number of iterates of (one forward and one backward). 

3. Finite representable multivalued maps 

For the proof it is necessary to use a computer (which can manipulate only a finite 
set of objects) to obtain cohomological information generated by the continuous 
map g . Since g is Lipschitz, the simplicial approximation theorem guarantees 
that in principle this can be done but the computations must be performed in 
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such a way that they apply not only to the simplicial approximation but also to all 
nearby maps, including the original map g . Hence it is more natural and technically 
simpler to employ finite multivalued maps. The ideas behind this solution will now 
be explained (sec [6] for details). 

On a computer one can work with infinite sets of mathematical objects (numbers, 
vectors, functions, relations, etc.) only by choosing a finite subset of objects and 
coding its elements with natural numbers. The elements of the selected finite subset 
of objects are referred to as representable objects (relative to the selected coding). 
In the proof of Theorem 1 the set of representable real numbers is given by the 
standard floating point representable coordinates. The representable sets in R'^ 
were chosen by selecting a compact set M C within which the dynamics of 
interest occur, a set D of representable vectors, and a representable number 77 , in 
such a way that the set M is covered by the collection of balls B{d , rj) , d E D . 
Since the particular metric chosen influences only the efficiency of the algorithm, 
for the proof the sup norm was chosen. Thus, the balls are in fact cubes. A 
representable set is any union of a subcoUection of these cubes. Let Mq C M . A 
representable multivalued map on Mq is a multivalued map : Mq V{M) such 
that the set {!F{x) \ x S Mq\ is a finite collection of representable sets. 

Representable multivalued maps allow for the passage from continuous dynamics 
of the Lorenz equations to the discrete dynamics of the computer. Furthermore, 
they are more readily attainable than simplicial approximations, as can be seen 
from the following description. Let / : A^q ^ M be a Lipschitz continuous function 
with Lipschitz constant L . Assume that for every representable vector d one can 
compute a representable vector /o(rf) such that ||/(rf) — /o(c^)|| < 5 , where 5 is a 
given representable number. Defining !F{x) as the smallest convex representable 
set which contains i?(/o(a;) , d + Lrj) , one obtains a multivalued map ^ : Mq M 
such that f{x) G J^{d) for all x G B{d , ??) . Then 

J^^'ix)- U :F(d), and jP'(a;) := f| JF(d) 

||d-x||<77 l|d-x||<77 

are easily seen to be finite, representable, and respectively upper and lower semi- 
continuous extensions of / . Furthermore, J^^ is convex valued. 
Letting 6^0 results in the following proposition. 

Proposition 5. Let Mq C R" be a compact set, and let f : Mq — > R" be a Lips- 
chitz continuous function. Then there exists a sequence J^n of finite representable 
extensions of f such that ^ f as n 00 . 

4. Numerical computations 

Simple and standard numerical integration of the Lorenz equations at the param- 
eter values (45, 54, 10) leads to a return map which strongly suggests the existence 
of horseshoe dynamics. In particular, one can find two rectangles Rq and i?i in 
the plane P which, under the numerical integration, appear to cross themselves 
and each other transversally. Thus, one expects the existence of an invariant set 
which is conjugate to the full shift dynamics on two symbols. 

As was suggested in the introduction, the proof of the semi-conjugacy described 
in Theorem 1 was obtained by computing a representable multivalued admissible 
extension Q of the Poincare map g such that N := Nq U Ni is an isolating 
neighborhood under Q . The sets A^o and A^i were taken to be carefully chosen 
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rectangles within Rq and Ri in order to minimize computation. Ideally the multi- 
valued extension Q would have been obtained by integrating the Lorenz equations 
numerically from the center of each cube in the grid of representable sets and in- 
corporating all errors into the size of the assigned value. Unfortunately the growth 
of errors is exponential with time, which caused technical difficulties that needed 
to be overcome. 

There are four sources of error: 

(i) the approximate arithmetic of the machine; 

(ii) the numerical procedure used to integrate the equations; 

(iii) extending the value from the center of a cube to the whole cube by means 
of Lipschitz constant estimates, as was explained in the previous section; 

(iv) estimating the point of intersection of the trajectory with the cross-section 
from two consecutive steps of the numerical method. 

In every case rigorous error bounds can be obtained (see [4] for details). It should 
be mentioned here that double precision arithmetic and the standard fourth-order 
Runge-Kutta method with step size 100/2^'^ were used. The error estimates for 
(iii) were based on the Gronwall inequality, the local Lipschitz constants, and log- 
arithmic norms. The second-order Taylor expansion of the solution was used to 
estimate errors in (iv). This resulted in the growth factor (ratio of the size of a 
value to the size of the grid) for equation (1) and in the function Q described 
above being approximately 10^ . The number of initial grid points needed to fulfill 
condition (2) by the associated multivalued map Q is proportional to the square 
of the growth factor. 

With the computing power available a significantly smaller growth factor was 
needed. To achieve this, twenty-three intermediate cross-sections, labelled Ek , k = 
— 1 , . . . , 22 , were introduced, with the original Poincare section appearing as P = 
E!_i = S22 ■ Using the technique described above, twenty-two finite representable 
multivalued maps Qi : S^^i — ^ were obtained. Observe that each multivalued 
map is an extension of the flow deflned map gi : S^-i —>■ determined by the 
Lorenz equations. Q was taken to be the composition Q21 o ^20 o . . . o Qq . 

The calculation was performed beginning with approximately 700,000 cubes cov- 
ering A^i . (It was not necessary to compute Q on iVo because of the symmetry 
present in the Lorenz equations.) The growth factor in this case was approximately 
30, which resulted in the diameter of Q{q) being less than 0.044 for all q ^ Q . 
Thus to determine whether N is an isolating block under Q , one need only check 
that 

B(g*-^(7V) n TV ng(iV), 0.044) cint N. 

This is indeed the case. Therefore N is an isolating neighborhood of Q , and hence 
by Theorem 3 is an isolating neighborhood of g . 

5. The Conley index of Q 

Recall that to compute the Conley index of an isolated invariant set for a map, 
one first needs to find an index pair (A'^ , L) and then to determine the index map, 
i.e. 

75 : H*{N, L) H*{N, L). 

For this particular index pair these objects can be determined essentially by in- 
spection, since the only cohomology group of interest is H^{N , L) . However, for 
more general problems it should be observed that D can be viewed as the set of 
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vertices of a simplicial decomposition of P , and hence Q is a multivalued simpli- 
cial map. Therefore, determining Ig is in fact a finite computation and hence in 
principle can be done by the computer. Prom Ig one computes the Conley index 
Con*(Inv(A^ , Q)) and hence by Theorem 3 the index of interest Con*(Inv(iV , g)) . 
For this problem the index computations lead to the algebraic hypotheses of The- 
orem 2, and hence Theorem 1 follows. 
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